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Abstract 

The non-Markovian nature of polymer motions is accounted for in folding ki- 
netics, using frequency-dependent friction. Folding, like many other problems 
in the physics of disordered systems, involves barrier crossing on a correlated 
energy landscape. A variational transition state theory (VTST) that reduces 
to the usual Bryngelson- Wolynes Kramers approach when the non-Markovian 
aspects are neglected is used to obtain the rate, without making any assump- 
tions regarding the size of the barrier, or the memory time of the friction. 
The transformation to collective variables dependent on the dynamics of the 
system allows the theory to address the controversial issue of what are "good" 
reaction coordinates for folding. 
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According to the energy landscape theory [0 , protein folding can be seen as a stochastic 
motion of a few collective coordinates describing protein conformation, on an average ther- 
modynamic potential To first approximation this motion is Brownian, and the folding 
time can be computed from diffusive rate theory A good quantitative comparison be- 
tween the analytical energy landscape theory and lattice simulations of smaller proteins has 
been made 0]; The diffusive behavior of the reaction coordinate's motion is approximate; 
in lattice simulations a fraction of the trajectories have ballistic crossings over the barrier, 
while others are quite diffusive suggesting a wide range of time scales for the collective 
reaction coordinate. Non-Markovian dynamics, expected on a rugged energy landscape, will 
affect reaction rates when the time to cross the top of the barrier is comparable to the 
memory time of the fluctuating forces acting on the collective coordinate. Analysis of such 
situations for reactions in condensed phases has led to a number of good approximations 
for rates 0-0. In our treatment we use variational transition-state theory (VTST) [§-10 



to discuss how non-Markovian dynamics of the chain affects folding, and to address the 
question of what is the best reaction coordinate for folding. We apply our results to the 
motion of a heteropolymeric protein chain, but many of the same issues occur for kinetics 
in other disordered systems, for example the nucleation of a crystal from a glassy liquid. 
We first discuss the form of the effective frequency-dependent friction ({w) for motion on a 
correlated energy landscape, specifically for a heteropolymer. Then we apply VTST to find 
corrections to the Kramers folding rate due to memory effects in C{uj), and anharmonicities 
in the potential. 

Protein conformational motion can be mapped onto a generalized master equation, with 
escape rates from a statistical ensemble of configurations given in terms of a waiting time 
distribution Pq{t, T) The configurational states may be grouped together in strata with 
a common value of their similarity to the native conformation, which is an approximate 
reaction coordinate for folding. This is often taken to be Q |]TIi|5|, the fraction of native 
contacts, but other choices are possible. Finding a "best" reaction coordinate is currently 
of great interest [p|,p]2|-[T^ . By projecting onto this coordinate, a diffusion equation for 
its probability distribution is obtained with a frequency- dependent diffusion coefficient [|l|. 
Correspondingly, the coordinate's motion can be characterized by an overdamped generalized 
Langevin equation (GLE): —dF{Q)/dQ — jQ (ir^Q(t — r, T)Q{t) +^(t) = , with a frequency- 
dependent friction coefficient (q{uj,T) = k^T j D^^iuj) satisfying {i(t)^(t')) = T(Q(t — t',T). 

The GLE implies that Q responds linearly to fluctuations in the other coordinates of 
the polymer chain apart from the nonlinearity inherent in the thermodynamic potential 
for Q. This should be a good approximation if many individual configurational states of 
the polymer chain are sampled for each value of Q, as expected above the glass transi- 
tion temperature of the stratum at Q. (q{u,T) is given by averaging over Pq{t,T) as 
Cq(cu,T) = Ag(r/(l + cur))/(l/(l + cur)) or / (^t^{^~''^) where A is the Laplace 

transform. The conformational motion distance scale is set by Aq = 2k^T / [AQ'^'-Jq) , where 
/S.Q is the step size and 7q is the probability a jump changes Q. Pq{t, T) on a correlated en- 
ergy landscape is obtained by first defining a local progress coordinate q as similarity to the 
given trap state, so the typical rate of escape |jl5| involves the calculation of the free energy 
barrier (FqIq^) — Fq(1)) for motion to states having native similarity near to Q, assuming 
the elementary moves are sufficiently local. A bilinear approximation to the entropic part of 
FqIq) can be used |15|, since contacts formed at small q, i.e. for a more weakly constrained 
polymer, cost more entropy than for a strongly constrained one at high q. The escape time 
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oc exp(AF^/T) depends on 2 parameters: 1.) The reduced temperature T = T /Tg-, where 
Tg = {AEq/2SqY^'^ is the glass temperature for the Q stratum. Sq is the configurational 
entropy at Q while AEq = (1 — Q){AEIj + QAE%) is the energetic variance of the states 
in terms of the variance of native (N) and non-native (M) contacts, and 2.) The reduced 
energy of the trapped state E = E/Egs, where Egs = —{'^SqAEqY^'^ is the ground state 
energy of the ensemble of states at Q. The escape time from a trap with energy E, '^q{E, T), 
is given by roexp[S'q(l — q^)(2E/T — \/a^ — 1/T^)6'r]. Here is the location of the barrier 
peak, and a"'" = (1 — q^)SQ/S'^, where S'^ is the fraction of Sq at the barrier peak (for the 
64-mer gt ^ 0.3 and ^ 1.6, see fig. 4-9 of ^). For states with E < EA{f) or when 
T > Ta{E) (determined by setting the barrier to zero), escape becomes downhill with short 
life-times of roughly the Rouse-Zimm time(s) ~ r^, hence 6r = 6{E — Ea{T))6{Ta{E) —T) in 
the exponent. At T = 1 where the system is frozen into one of a few ground states {E = 1), 
and the typical escape time is tq{1,1) = Toexp[SQ{l — q^){l — — Tq exp(0.27 S'q^^'*). 

Correlations lower the free energy barrier to correspond to roughly 1/4 the total configu- 
rational entropy |]T5[, as opposed to the full entropy as in uncorrelated landscapes. The 
distribution of occupied state energies E at temperature T is a Boltzmann weighted gaus- 
sian: Pq{E,T) ~ exp[— S'q(-E — 1/T)^]. Reflecting this, the distribution of escape times is 
easily calculated as PQ(r,f) = j\dE PQ{E,f)6[T - TQ{E,f)] = P^"^ + P^. Apart from 
barrier-less escapes with P^~^{t,T) = A{T)6{t — Tq), this yields essentially a log-normal 
distribution: 
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where 9^ = 9{TA{l) — T)9{T—To)9{TQ{l,T)—r), S = Sq, and A = Aq(T) is the normalization 
constant of Pq{E,T), while A = A(r) is the weight for fast escapes from states without 
barriers. 

The above analysis applies at temperatures higher than the thermodynamic glass tran- 
sition temperature Tq (but below the activation temperature Ta, where escape from all 
occupied states is barrier-less). At or below Tq the temperature independent distribution 



of state energies becomes P{E) ~ exp 
distribution of escape times P{E{t))/ 



'E/Tg). Arrhenius-like escape from states yields the 



^ f-i^+T/Ta)^ giving stretched exponential relax- 
ations |]I6|,|17|. In this regime the assumed linearity of Q dynamics is questionable however. 

Evaluating Dq{u!,T) or (q{lj,T) 

above Tq, note that (e~*/^) ~ rfze"^^" exp[-6^e-"25^] (^here ca = 2(1 - q^)/T, and 

b = exp[S'(l — q^){l/a'^ — 1/T^)]) is reminiscent of the after-effect function f{bt/To,C2/VS) 
used for non-exponential decay in glasses [|TB|. But for the mesoscopic systems relevant to 
folding (A^ ^ 100), a better approximation is to linearize the gaussian term on the range 
(zl, Zu) = (T/2a^ — 1/2T, 1 — 1/T), yielding a closed form expression for the friction 

^ 1^ + ig^ ir-'J". (1, 1 - c, 2 - c, -cu/f) (1, 1 - c, 2 - c, -g;)] 
^ + f (1' -c, 1 - c, -^) - /^^ (1, -c, 1 - c, -u/f)] 



Cq(^,t) 



where 2^ is the hypergeometric function, u is in units of I/tq, F = Ae^^'^ / C2b^S , c = 
2zm/c2, f = To/tq{\,T) < 1, and = {z^ + 2l)/2. The corresponding diffusion coefficient 
Dq{uj,T) = k^T / Cq{ijj ,T) is plotted in fig. 1. 
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Non-Markovian rate theory proceeds by recognizing that the GLE is equivalent to a 
particle bilinearly coupled to a bath of oscillators |]19| with an effective Hamiltonian Ti. = 
pl/2 + F{Q) + lJ2j[pl^ + {i^jXj -CjQ/ujf] with ({t) = T,j{Cj/uJj) cos{ujt). These additional 
collective modes describe the dynamics of Q fluctuations within linear response. Rather than 
dealing with non-Markovian dynamics of Q we can study the dynamics of this equivalent 
many-dimensional system without memory. When a single barrier exists in F{Q) it makes 
sense to use multidimensional transition state theory. If the barrier is large and its parabolic 
part dominates, Ti is quadratic and may be diagonalized by a normal mode transformation, 
which singles out as a reaction coordinate an unstable mode p with imaginary frequency iX^ 
given by the solution of A%(3t(A^,T) = muj^'^ = \d'^F{Q)/dQ'^\QX. This frequency is identical 
to the (overdamped) Grote-Hynes reactive frequency 0]. The effect of friction is to leave the 
barrier height unchanged, but to rotate the reaction coordinate to a different direction in 
configuration space. When the friction has no memory and the Q motion is purely diffusion, 
the reactive frequency is that of an overdamped inverted harmonic oscillator corresponding 
to the Kramers prefactor in the rate. For a general F{Q) we can separate the quadratic 
part from the anharmonic part of F . Then in the normal coordinates at the saddle point 
H. = \[p'^p — A^^p^ + J2j{plj + ^'jV])] + Fi{uooP + Y.j UjoVj)- The coefficients Ujo are elements 
of the orthogonal normal-mode transformation such that q = Q — Q^ = UooP + Z)j '^joUj ( 
is given in terms of the friction kernel as ul^ = l/[CQt(A*,T)/A^ + d(Qt{s,T)/ ds\^^^x])- 
addition to the Grote-Hynes coordinate p one can define a residual collective bath coordinate 
(7 = (1 — Uoo)^^/^ Y.j UjoVj = Y.j UjoVj, which appears in the anharmonic part of the 

potential. The effects of dynamic friction, reflected by re-crossings in the Q coordinate, 
are in fact accounted for by ballistic motion across a new thermodynamically determined 
dividing surface involving a 2-D potential in (p, o") ||^. The equilibrium flux across any 
dividing surface written in these new coordinates is 

^_ Jdp, dp dp^ da 6 if) (V/ ■ p) e (V/ ■ p) e-^^' 
J dpp dp dpfj da e~^^* 

Here f = p — g{<y) serves as a new progress coordinate. / = determines a dividing surface 
between reactants and products, the factor 6{f) localizes the integration to that surface, 
while V/ ■ p = Ppdf/dp + p^df /da is the flux density across the surface. The reduced 
Hamiltonian depends only on two coordinates p and a: H} = -f- — A^^p^ + fl'^a'^] + 
FiiuooP + uia), with a collective bath frequency Vl? = u^/ {ul^/ X^"^ — l/cu^^). Carrying out 
the integrations, the rate can be written as a correction to the Grote-Hynes rate: F = 
P{Xyu!^)To. Tq = {uJo/2n) exp(— /^F-t) is the TST rate, while A^^^^ is the G-H transmission 
factor. The additional correction to the G-H prefactor is given by a quadrature: 




P = 

Here E[g] = ^{Q'^a'^ — X^'^g{ay)+Fi{uoog{cr)+Uia). If Fi vanishes, the potential is parabolic, 
and variational minimization of the rate yields / = p as the progress coordinate and q = 
Uia = a as the ideal dividing surface. The transition state position is coupled strongly to 
the bath mode, but the correction gives P = 1, reproducing the G-H rate F = {X'^ /uj^)Vq. 
For systems with large barriers ^ IQk^T , the parabolic approximation discussed above is 
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highly accurate but we can find corrections to it by finding a more general planar dividing 
surface / = UooP+J2j UojVj — po = (where po is the distance of the dividing surface from the 
barrier peak) oriented such that the TST rate is minimized. Using the variational framework 
developed in ^ along with the correlated landscape theory of (q{u,T), we find corrections 
to the GH rate at conditions of equilibrium between the folded and molten globule states (at 
Tp), see fig. 2. The optimal barrier location po (see inset A) is little different from the naive 
choice of the free energy maximum, for the 64-mer. The rate enhancement over the Kramers 
result has a peak at intermediate T^/Tq due to the larger relative frequency dispersion of 
the diffusion coefficient at intermediate temperatures. For systems this large, the planar 
dividing surface assumption is accurate to within 5%. 

For a 27-mer imitating a small protein at Tp = I.GTq, we can use the simulated auto- 
correlation function of Q, c{t) (see fig. 9 of 0]), to determine C('^)- C('^) = {Q^) / {1 /c{uj) — 
uj). Here c{uj) = Ctc{t). For the 27-mer at Tp, the potential is very anharmonic. A planar 
dividing surface is no longer optimal. For such low barriers, minimizing the fiux through 



the TST surface can be carried out using the calculus of variations as done by Miller |20 



and PoUak [|TU|. Finding the stationary point of the rate for arbitrary variations of the 
dividing surface functional Sg{a) yields a differential equation for g{(T): g"/{l + g'"^) = 
j3{g'dE[g]/ da — dE[g]/ dg). Treating g and a as parametrized variables in terms of an 
independent parameter t (such that g' = g/&, and setting the first integral |((?^ + cr^) = 
Eq — Vib) recasts this variational equation into Hamilton's equations of motion for g and a on 
an effective temperature dependent potential Vg = — {1/2 P) exp{— 2 PE[g]) at total energy 
|(Pg+p^) + Vg = 0. The optimal dividing surface is a classical periodic trajectory on Vg with 
infinite period, that divides the (p, a) space into reactants and products. The correction P 
is given in terms of the action along this optimal trajectory: P = (/3^r2^/27r)^/^ / ds ^J—2Vp 
where ds is arc-length along the trajectory. The optimal dividing surface for the 27-mer is 
plotted on the potential E{p, a) in fig. 3. The correction P = 0.85. The rate is moderately 
reduced from the G-H value of 1.57 k^^, giving k = 1.33 k^^ for the corrected rate, in closer 
agreement with the naive Kramers value. While the Kramers approximation made in 
is numerically quite accurate, the optimal dividing surface shows the transition region is 
quite spread out in the coordinate Q. By finding the optimal dividing surface, VTST seeks 
that coordinate which behaves most ballistically. Including the Q dependence of ( at higher 
nativeness may explain some of the discrepancy between simulations and rate calculations . 

The present analysis can easily be generalized to include ordering along additional col- 
lective coordinates, as for example the total density of contacts, which is often important 
if collapse is not fast |^. Potentially interesting effects may arise in such scenarios due 
to anisotropic friction |^, and are a topic of future work. The methods presented here are 
general, and also apply to other condensed matter systems with rugged landscapes, e.g. 
glasses and clusters. Nonlinear couplings between the system and bath may also be treated; 
This may allow explicit treatments of the deep traps in the glassy regime. 
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FIGURES 




CO To 

FIG. 1. ToDQt{uj,T) for 64-mer near the transition state w 0.3, as a function of frequency 
(in units of To), at several temperatures between the thermodynamic glass temperature and kinetic 
(activated) glass temperature. There is a rapid increase from the zero frequency value Dq{0,T) 
dependent on the typical escape time, to a higher asymptotic value depending on how many of the 
states are untrapped and have short lifetimes at that temperature. The dispersion in the values 
of the diffusion is thus maximum at intermediate values of temperature. The largest values of the 
diffusion constant are set by ctq = k^T/X^, and at T4, DQt{uj,TA) = (J^= AQ'^^q/2to ~ 0.0015/ro 
for the 64-mer. 
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FIG. 2. 

(solid line) Ratio of the Grote-Hynes rate to the Kramers rate, as a function of T^/Tq. 
(dashed line) VTST rate enhancement including anharmonic effects of a finite size barrier. 
The effect here is small, however larger effects are seen for shorter polymers (see text). The 
system is a 64-mer at folding equilibrium. The variance of interaction energies is varied, so 
that the temperature ratio T^/Tq varies {Ty/Tq = ^/rj + \/'r] — 1, with r] = El/2SoAEl.j, 
where E-^ is the native state energy, and So is total entropy). The rate closely follows the G- 
H result k'^^'/k™ = CQt(0, T)/CQt(A™, T). (Inset A) Thermodynamic potentials vs. Q. For 
more rugged landscapes the barrier is flatter, and this reduces the prefactor to the rate since 
there are more recrossings. The Tp/T^ values are 2.45, 1.84, and 1.24 in order of decreasing 
barrier size. The small vertical bars near the barrier peak are where the VTST dividing 
surfaces cross the coordinate Q. (Inset B) Enhancement of the rate at the maximum 
value of 3.8 at Tp/Tg ~ 1.6 by allowing the prefactor Aq = k^T/a^ in D(A°",T) to vary 
(ctq 0.0015 is the original value). 
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FIG. 3. (Color) Potential surface E{q,uia) (red) and in the normal coordinates {uooP,uia) 
(inset), along with the variational dividing surface which minimizes the TST flux (heavy line). 
Contours are drawn at intervals of about 2k^T. The potential in the Markovian case (dashed 
green), with the corresponding Kramers rate, is further skewed with respect to the dividing surface, 
indicating paths in this case are even more diffusive. X marks the position of the molten globule 
minimum, and O the native minimum. The short vertical lines (and horizontal lines in inset) 
bound a region of ^ 70% of the total flux across the dividing surface. In {q,uia) space there is 
flux over a wide range of q values, A.q/{qnative — Qmg) ~ 0.44, so that the transition state theory 
that reproduces the multiple crossing physics in Kramers theory does not have a well defined value 
of g-f. However in {uooP,uia) space the TST surface tends towards orthogonality to the reaction 
coordinate UqoP'- ^p/ {pnative — Pmg) ~ 0.04, indicating trajectories behave more ballistically along 
P- 
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